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§1.  OBJECTIVES 

The  primary  goal  of  the  project  was  the  development  and  proof  of  concept  of 
a  biologically  inspired  topology  optimization  method  for  aerospace  design. 

§2.  ACCOMPLISHMENTS 

The  objective  set  forth  for  this  project  has  been  fully  attained.  Indeed,  the 
biologically  inspired  design  methodology  developed  in  this  project,  hereafter  named 
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bioTOM,  has  been  thoroughly  assessed  through  a  large  number  of  test  cases.  The 
test  cases  ranged  from  classical  benchmark  problems  in  topology  optimization  to 
realistic  engineering  problems.  These  test  cases  proved  the  concept  and  showed 
the  well-suitability  of  bioTOM  for  aerospace  design. 

Next  in  this  section,  a  brief  account  of  the  developed  design  methodology  is 
presented.  This  is  followed  by  a  summary  of  the  test  cases  studied  in  this  project. 
The  section  ends  with  the  main  findings  of  the  research. 

1.  bioTOM 

As  Richard  Dawkins  pithily  puts  it  “Life  is  the  execution  of  programs  written 
using  a  small  digital  alphabet  in  a  single,  universal  machine  language”  [1],  The 
programs  that  indirectly  determines  the  entire  developmental  process  of  living 
organisms  are  encapsulated  in  our  DNAs.  And  the  latter  are  evolved  according  to 
Darwin’s  laws  of  natural  selection.  Exploiting  this  genetic-evolutionary  model  for 
engineering  is  the  novel  paradigm  set  forth  in  this  work. 

This  sub-section  starts  with  an  outline  of  the  elements  of  map  L  systems  that 
are  relevant  to  the  present  work.  The  reader  interested  in  the  general  theory  of 
map  L  systems  may  refer  to  [2]  and  the  references  in  that  book. 

At  the  end  of  f960’s,  the  eminent  biologist  Aristid  Lindenmayer  introduced 
a  novel  type  of  grammar  system,  where  the  rewriting  was  carried  out  in  parallel 
[3,  4],  rather  than  the  sequential  substitution  in  previous  formal  languages  [5].  The 
parallel  rewriting  system,  thereafter  named  Lindemayer  system  (or  L  system  for 
short),  was  initially  proposed  to  model  the  development  of  the  branched  topology 
in  plants.  Over  the  years,  however,  L  systems  became  an  important  scientific 
theory  in  itself  [2] ,  with  multiple  ramifications  ranging  from  the  original  model  for 
plant  development,  to  the  mathematics  of  formal  languages,  to  computer  graphics, 
to  music.  Map  L  systems  [2,  6,  7]  extend  the  parallel  rewriting  in  L  systems  to 
planar  graphs  with  cycles,  called  maps  [8].  The  maps  are  developed  according 
to  cellular  division  rules  that  have  originally  been  devised  to  model  the  cellular 
division  process  in  simple  organisms. 

Formally,  a  map  is  defined  as  a  finite  set  of  regions.  Each  region  is  bounded 
by  a  sequence  of  edges  and  the  edges  intersect  at  vertices.  Every  edge  is  part  of 
the  boundary  of  a  region  and  the  regions  are  simply  connected.  These  maps  are 
analogous  to  cellular  layers,  where  the  regions  represent  the  cells  and  the  edges 
their  walls.  To  develop  these  maps,  the  Binary  Propagating  Map  OL-systems  with 
markers  or  mBPMOL-systems  proposed  by  Nakamura,  Lindenmayer  and  Aizawa 
[2,  7]  is  selected — mBPMOL-systems  are  more  powerful  than  similar  context-free 
map  L  systems  [2],  The  method  is  called  binary  because  during  the  cell  division 
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the  cells  divide  in  two.  It  is  propagating  since  cells  cannot  fuse  or  vanish.  The 
designation  OL  systems  refers  to  context-free  parallel  rewriting  systems  that  do 
not  allow  for  region  interactions.  Finally,  markers  specify  juncture  points  at  the 
edges  where  the  cell  can  divide — the  markers  are  analogous  to  attachment  sites 
for  division  walls  during  mitosis  [2]. 

Mathematically,  mBPMOL-systems  consist  of  an  alphabet  E,  an  axiom  u>  and 
a  finite  set  of  rules  P.  The  alphabet  is  a  finite,  non-empty  set  E,  whose  elements 
are  called  letters  or  tokens.  An  alphabet  can,  in  general,  contain  any  symbol  as 
letters.  A  simple  example  of  an  alphabet  would  be  E  =  {A,  B,  C, ...,[,  ],  +,  — }.  An 
axiom  is  a  non-empty  word  of  letters  from  the  alphabet.  The  number  of  letters  in 
the  word  should  be  equal  to  the  number  of  edges  in  the  initial  map.  For  example, 
using  the  previous  alphabet  and  for  an  initial  map  of  four  edges,  u  =  AABC 
is  an  axiom.  Each  rule  is  of  the  form  A  — >  a,  where  the  edge  A  £  E  is  called 
the  predecessor,  and  the  string  a,  composed  of  symbols  from  E  including  special 
symbols  [,  ],  +  and  — ,  is  called  the  successor.  Using  again  the  previous  alphabet, 
an  instance  of  a  rule  would  be  B  — >  A[—C}.  So,  when  rewriting  a  word,  every 
occurrence  of  the  letter  B  would  be  replaced  by  A[—C\. 

The  symbols  [  and  ]  specify  makers  for  possible  cell-dividing  walls.  Symbols 
outside  of  the  square  brackets  specify  the  number  and  type  of  edge  subdivisions — 
each  subdivision  with  the  same  length.  Inside  the  brackets  the  first  symbol  is 
either  +  or  — ,  placing  the  marker  to  the  left  or  to  the  right  of  the  predecessor 
edge,  respectively.  The  second  symbol  within  brackets  is  always  a  letter.  The 
letters  can  carry  an  arrow  over  them  to  represent  the  local  edge  orientation  of 
the  successor  edges  relative  to  the  predecessor  edge.  An  example  of  a  rule  with 

orientation  is:  A  — >  BA[—C]x. 

The  cell  division  process  is  effected  in  the  derivation  phase,  where  each  cell  is 
scanned  searching  for  matching  markers.  If  in  a  cell  there  exist  two  markers  in 
different  edges,  and  both  markers  carry  the  same  letter  and  follow  along  the  same 
orientation,  then,  depending  on  division  criterions  explained  below,  a  cell  division 
can  be  formed  connecting  these  two  markers.  More  than  one  pair  of  matching 
markers  can  be  found  in  one  cell,  however.  In  this  case,  and  owing  to  the  binary 
character  of  the  method,  only  one  pair  of  markers  is  selected  for  connection  and 
the  remaining  markers  are  discarded. 

The  selection  and  validation  of  the  edge  pair  is  performed  according  to  two 
cell  division  criterions:  (1)  the  angles  in  the  divided  cells  must  be  larger  than  a 
prescribed  lower  limit;  (2)  the  areas  of  the  offspring  cells  must  be  larger  than  a 
predefined  percentage  of  the  parent  cell.  The  first  criterion  prevents  the  creation  of 
cells  with  too  narrow  wedges  angles,  while  the  second  precludes  excessively  slender 
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or  relatively  minute  offspring  cells  to  form.  When  multiple  markers  are  available, 
the  first  pair  to  pass  the  two  criterions  above  are  selected  and  the  rest  are  dropped. 

The  cellular  division  process  explained  above  is  best  grasped  with  an  example. 
Consider,  for  instance,  the  alphabet  E  =  {A,  B,  x,  [,],+,  — },  the  initial  map  or 
axiom  uj  =  ABAB  and  the  production  rules: 

A  -»•  B[-A]x[+A]B 
B  -*■  A 


The  symbols  for  which  a  rule  is  assigned — A  and  B  in  this  example — are  called 
non-terminal  tokens,  whereas  the  remaining  symbols  are  called  terminal  tokens. 
The  rules  for  these  latter  tokens  are  not  listed  because  these  tokens  are  constants; 
for  instance,  a  symbol  x  in  a  word  is  copied  as  x  in  the  rewriting. 

The  derivation  process  starts  with  each  edge  of  the  initial  map  being  assigned 
a  label — see  figure  1  top — and  a  global  orientation  that  defines  the  left  (+)  and 
right  (-)  branching — the  orientation  is  counter-clockwise  in  this  example.  Then 
each  edge  is  relabeled  and  divided  according  to  the  corresponding  rule.  Thus,  in 
this  example,  the  edges  labeled  B  are  relabeled  A  and  no  sub-division  takes  place. 
The  edges  labeled  A  are  sub-divided  into  three  new  and  equal  segments.  Consider, 
for  example,  the  edge  A  at  the  bottom.  The  first  segment  is  labeled  B ,  the  second 
x  and  the  third  again  B.  After  the  first  segment,  a  marker  A  ought  to  be  attached 
to  the  right  (the  symbol  —  precedes  the  marker  A  within  the  first  parenthesis), 
but  since  the  marker  would  lay  outside  of  the  map  this  action  is  voided.  Moving 
along  the  edge  according  to  the  global  orientation,  at  end  of  the  second  edge  x,  a 
marker  A  is  placed  at  the  left  of  the  edge.  Applying  this  rewriting  process  for  the 
remaining  A,  two  markers  are  obtained,  labeled  A  pointing  to  inside  the  initial 
map.  Both  markers  are  of  the  same  type  and  assuming  the  two  aforementioned 
criterions  are  satisfied  the  cell  is  divided  by  connecting  the  attachment  sites  of  each 
edge.  The  result  is  the  new  map  in  figure  1 — first  step.  Iterating  this  cell  division 
process  generates  a  sequence  of  maps  that  models  the  developmental  stages  of  the 
structural  topology. 

Remark  1.  Note  that  the  cell  division  would  stop  if,  at  a  certain  stage,  all  edges 
were  labeled  with  a  terminal  token.  However,  stopping  only  when  all  letters  are 
terminal  tokens  could  lead  to  infinite  iterations:  for  example,  if  x  would  be  absent 
from  the  previous  rules.  Thus  a  maximum  number  of  cell  divisions  is  specified  and 
the  actual  number  of  developmental  stages  for  the  topology  generation  is  optimized 
in  the  evolutionary  algorithm. 
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Remark  2.  The  cellular  division  process  described  above  corresponds  to  the  orig¬ 
inal  map  L  systems  by  Nakamura  and  Lindenmayer  and  Aizawa  [2,  6,  7].  Here 
it  is  modified  slightly  to  allow  for  variations  of  properties,  such  as  edge  radius, 
within  the  topology.  The  approach  is  analogous  to  changing  the  state  in  a  tur¬ 
tle  interpretation  of  L  systems  [2]  and  essentially  it  provides  a  way  to  modify 
the  characteristics  of  the  structure  with  rules.  So,  terminal  tokens  are  introduced 
that,  when  present  in  a  rule,  will  change,  by  a  quantum  amount,  some  specific 
properties  of  the  topology.  For  example,  the  symbols  *  and  /  specify  the  change 
of  the  radius  of  the  edge  by  a  quantum  ratio  Sr.  So  the  first  rule  in  the  previ¬ 
ous  example  could  now  include  these  symbols  as:  A  — >  *B[—  *  A\/x[+  *  A]B. 
With  this  new  rule,  at  the  second  developmental  stage  would  lead  to  the  map: 
ui 2  =  *A[—  *  *B[—  *  A]/:r[+  *  .4 1 B  /.r  [•  i  *  *B[—  *  A]/:r[+  *  A].E>]A  *  B[—  *  A]/:r[+  * 
A]B  *  A[*  *  B[—  *  A]/:r[+  *  A]R]/:r[+  *  *B[—  *  A]/:r[+  *  A]F>]  A  *  B\—  *  A]/:r  [+  *  A]B 
and  the  first  edge  A  would  have  a  radius  5r  *  r0  where  r0  is  an  initial  radius,  the 
first  edge  B  would  have  a  radius  *  ro  and  so  on.  In  this  way,  it  is  allowed  the 
control  of  the  edge  radius  along  with  the  developmental  stages.  Introducing  more 
tokens  of  this  type  allows  the  control  of  any  desired  design  variable. 

Remark  3.  As  explained  in  [2],  the  cellular  division  in  living  organisms  entail  two 
stages:  a  division  stage  proper  and  a  cellular  dynamic  stage.  So  after  each  cellular 
division,  actual  cells  deform  according  to  forces  acting  on  their  walls.  Because 
structural  topology  is  the  only  interest,  this  dynamic  stage  is  modeled  as  Picard 
iterations  for  the  equilibrium  of  two  forces  acting  on  the  map  vertices:  an  elastic, 
spring  force  on  each  edge  of  the  topology  and  an  “osmotic”  pressure  proportional  to 
the  inverse  of  the  area  of  each  region.  The  elasticity  constant  of  each  edge  and  the 
osmotic  pressure  constant  are  modified,  and  can,  therefore,  be  evolved  within  the 
topology,  using  the  state  rules  explained  in  the  previous  remark.  Also,  to  increase 
the  efficiency  of  the  overall  methodology,  the  maximum  number  of  Picard  iterations 
is  bounded  to  a  small  number. 

Remark  4.  The  original  map  L  systems  deals  with  convex  regions  only  [2,  6,  7].  In 
this  work,  topology  optimization  is  performed  with  regions  that  can  be  non- convex. 
So  to  extend  map  L  systems  for  maps  with  non- convex  regions,  a  third  criterion 
is  added  for  cellular  division,  where  it  is  required  the  new  edge  dividing  the  parent 
cell  to  remain  inside  the  parent  cell. 

At  the  end  of  the  division  stages,  a  planar  topology  is  defined.  With  this 
topology,  a  finite  element  (FE)  model  can  be  formulated  to  calculate  the  fitness  of 
each  design — the  details  of  the  FE  model  for  each  test  case  are  explained  below. 
The  fitness  of  each  design  can  then  be  evolved  using  an  evolutionary  algorithm. 


5 


Evolutionary  algorithms  [9,  10]  are  biological  metaphors  that  produce  high 
quality  designs  by  identifying,  recombining  and  enhancing  the  best  features  present 
in  an  adapting  pool  of  replicators.  Before  the  genetic  encoding  is  explained,  recall 
the  basic  elements  of  evolutionary  algorithms.  In  these  algorithms,  the  evolution 
starts  with  a  population  of  individuals,  each  of  which  carrying  a  genotypic  and 
a  phenotypic  content.  The  genotype  encodes  the  primitive  parameters  that  indi¬ 
rectly  determine  an  individual  layout  in  the  population.  Here  they  consist  of  the 
parameters  defining  the  map  L  system  and  the  geometry  and  physical  properties 
of  the  object.  The  phenotype  relates  to  the  ensuing  structural  model:  its  topology, 
geometry  and  physical  properties.  With  the  phenotype,  the  evaluation  of  the  fit¬ 
ness  of  each  individual  design  is  obtained  using  a  physico-mathematcal  model  and 
the  finite  element  method.  The  fitness  will  vary  from  problem  to  problem,  and  in 
the  case  of  multi-objective  optimization  multiple  functions  are  evaluated  for  each 
individual.  The  evolutionary  algorithm  then  advances  the  current  population  to 
the  next  generation  by  applying  selection,  mutation  and  crossover  operators. 

Genome  definition.  Fixing  an  alphabet  E,  there  are  three  classes  of  parame¬ 
ters  that  affect  the  topology  defined  by  the  map  L  system  in  our  variant  of  the 
mBPMOL-system.  The  first  is  the  axiom  word  u,  which  has  an  effect  on  the 
growth  and  dynamics  of  the  developmental  stages.  The  axiom  is  encoded  as  a 
word  defined  by  the  same  letters  of  E.  The  second  class  are  the  edge  production 
rules  P.  In  keeping  with  the  analogy  with  biological  systems,  the  production  rules 
may  be  seen  as  the  regulators  that  control  the  complex  processes  involved  in  the 
production  of  cells,  organs,  limbs,  etc.  by  parallel  interpreting  the  DNA.  Each 
production  rule  of  the  map  L  system  is  encoded  according  to  a  master  rule  of  the 
form: 

Y  — >  XxX2  . . .  Xn_iXn 

where  Y  is  a  non-teminal  token  and  A",  denotes  terminal,  non-terminal  or  special 
tokens.  The  number  of  slots  in  the  rules,  n,  is  decided  by  the  user  and  is  the  same 
for  all  rules — to  avoid  bracket  mismatches,  if  Xt  includes  a  bracket  X,  =  [±A,] 
is  defined  with  A,  a  terminal  token  or  the  non-terminal  token  x.  Finally,  the 
remaining  genome  is  filled  in  with  parameters  defining  the  geometry  or  physical 
properties  of  the  object  such  as  initial  edge  thickness  or  initial  edge  spring  stiffness 
modulus. 

In  summary,  the  genome  can  be  seem  as  a  partitioned  ribbon  whose  first  part 
is  occupied  by  the  letters  of  axiom,  the  second  by  the  production  rules  and  the 
third  by  geometric  and  physical  parameters  pertinent  to  the  body. 

Multi-objective  optimization  In  multi-objective  optimization  problems,  there 


6 


is  a  considerable  number  of  candidate  metrics  for  defining  optimality.  In  the 
present  work,  the  most  commonly  adopted  criterion  of  Pareto  optimum  is  chosen 
(see,  for  instance,  [11])  as  explained  below. 

The  optimization  problem  can  be  stated  as  follows:  Find  the  minimum  of 


subject  to 


g(%)  <  o, 


where  x  is  the  vector  of  design  variables,  fL  is  the  i-th  target  or  objective  func¬ 
tion  and  g  is  the  constraint  vector.  In  the  present  work,  constraints  are  enforced 
through  penalization  of  the  fitness  values.  A  vector  satisfying  the  constraint  is 
called  feasible.  A  feasible  vector  £  is  Pareto  optimal  or  nondominated  if  there  is 
no  feasible  vector  x  in  a  neighborhood  of  £  such 


fi(x)  </*(£) 


for  all  i  =  1, . . . ,  n  and  for  at  least  one  i  e  {1, . . . ,  n} 


fi(x)  <  /*(£)■ 


Hence  a  Pareto  optimum  is  a  point  where  around  it  you  cannot  measurably  improve 
some  targets,  without  simultaneously  worsening  others.  The  set  of  nondominated 
points  is  called  the  Pareto  front. 

With  the  criterion  for  optimality,  the  crowding  distance  is  chosen  to  compute 
the  fitness  function  [12].  The  crowding  distance  can  maintain  diversity  in  the 
population  and  has  good  convergence  properties. 

Following  the  description  of  the  biologically  inspired  method,  an  account  of  the 
test  cases  carred-out  in  this  project  is  presented  next. 

2.  Test  cases 

To  assess  the  suitability  and  efficiency  of  bioTOM  for  aerospace  design,  a  num¬ 
ber  of  test  cases  were  selected  and  simulated.  The  results  obtained  in  this  assess¬ 
ment  are  summarized  in  the  following  sub-sections. 

2.1.  Cantilever  benchmark 

The  short  cantilever  problem  is  a  classic  benchmark  for  structural  topology 
optimization  [13,  14,  15].  In  this  problem,  a  2  x  1  cantilever  in  plane-stress  and  with 
a  load  applied  at  the  free  extremity  (figure  2)  is  optimized  for  minimum  weight. 
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A 


Figure  1:  First  four  steps  in  the  cellular  division  process  modeling  the  develop¬ 
mental  stages  of  the  structure  topology. 

The  Young  modulus,  E ,  the  point  load,  P,  the  density,  p  and  the  cantilever  width 
w  are  all  set  to  lm,  while  the  Poisson  ratio,  u,  is  set  to  0.3.  These  values  do  not 
intend  to  represent  any  specific  material,  but  rather  to  provide  a  simple,  yet  well 
defined,  benchmark  problem  to  compare  the  performance  of  new  methodologies 
against  other  established  approaches. 

In  words,  the  objective  of  this  benchmark  is  to  obtain  the  lightest  cantilever’s 
structure  for  a  prescribed  maximum  deformation.  Mathematically  the  problem 
can  be  expressed  as: 


minimize  IT(x)  (1) 

such  that  Dmax(x.)  -  DUm  <  0 

where  IF(x)  is  the  cantilever’s  weight  nondimensionalized  by  the  weight  of  the 
rectangular  cantilever  2x1.  Dmax(x )  is  the  maximum  displacement  and  Ditrn  is 
maximum  displacement  allowed,  in  this  case  Dum  is  set  to  220.  Finally  x  is  the 
vector  that  encodes  the  Map  L-system  topology  as  explained  previously. 


Figure  2:  The  2x1  cantilever  benchmark. 


The  fitness  value  is  evaluated  with  plane  stress  finite  elements  available  in 
the  FE  toolbox  COMSOL.  This  software  supports  a  scripting  language  similar  to 
MATLAB  that  is  very  useful  to  automatize  the  calculation  of  the  fitness.  Geometry 
creation,  meshing,  solving  and  postprocessing,  all  can  be  done  through  scripting. 
Some  of  the  commands  used  are  listed  bellow. 

•  curve2,  creates  a  2D  curve  object  in  the  form  of  a  Bezier  curve; 

•  elipt2,  creates  a  solid  ellipse/circle; 

•  meshinit,  meshes  the  geometry  using  the  Delaunay  triangulation; 

•  femstatic,  solves  the  stationary  PDE  problem  with  a  nonlinear  or  linear 
solver; 

•  postint,  integrates  expressions  over  subdomains,  boundaries,  edges,  and 
vertices; 

•  postmin,  computes  the  minimum  value  of  an  expression. 

A  complete  reference  for  the  available  commands  can  be  found  in  the  COMSOL 
user’s  manual. 

This  benchmark  problem  has  one  constraint  and  as  mentioned  previously  con¬ 
straints  are  enforced  through  penalization  of  the  fitness  values.  Thus,  the  penal¬ 
ization  Pw  is  added  to  the  fitness  value  in  equation  2.  Pw  is  given  by: 


p ^  _  (J  ^pnax(Drnax  (x) — Dj jm  ,0)  _ 


(2) 


with  (7  =  1; 
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The  convergence  history  for  the  best  individual  in  a  generation  is  shown  in 
figure  3.  Initially,  the  improvement  between  generations  is  fast  but  eventually  it 
flattens  out  as  the  best  designs  become  more  refined.  After  60  generations  the 


Figure  3:  Convergence  of  the  bioTOM  for  the  minimum  weight  optimization. 

lightest  cantilever  geometry  that  satisfies  the  optimization  2  weights  W  =  0.26. 
This  means  that  with  a  proper  structure,  a  cantilever  with  26%  of  the  mass  of  the 
rectangular  cantilever  from  figure  2  exhibits  a  maximum  deformation  smaller  than 
Diim.  The  structure  and  the  Von  Mises  stresses  for  this  design  are  shown  in  figure 
4.  From  these  figures  it  can  be  seen  how  the  material  is  placed  such  that  there 
are  no  unnecessary  elements  in  structure,  i.e.,  all  elements  withstand  similar  level 
of  stress.  Also,  it  can  be  seen  that  the  elements  that  are  clamped  are  the  thickest 
ones.  This  is  in  agreement  with  what  was  expected  given  that  at  the  cantilever’s 
root  the  stresses  are  maximum. 

Table  2.1  compares  the  results  obtained  with  bioTOM  against  other  the  results 
of  methods  found  in  the  literature:  using  a  bit-array  representation  [14,  15]  or  a 
Voronoi-based  representation  [13],  for  the  structure.  The  column  Wmin  demon¬ 
strates  that  the  optimal  structure  for  the  bioTOM  represents  a  significant  im¬ 
provement  with  respect  to  the  other  methodologies.  The  column  feasible  shows 
the  percentage  of  feasible  individuals  for  all  generations.  As  the  results  of  the 
table  shows,  bioTOMhas  a  vital  advantage  over  the  Genetic  Algorithmof  creating 
feasible  desings.  This  is  relevant  because  the  efficiency  of  the  genetic  algorithm 
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Figure  4:  (a)  The  optimal  cantilever’s  structure  for  the  minimum  weight  optimiza¬ 
tion  after  60  generations,  (b)  The  colormap  shows  the  Von  Mises  stresses  from  0 
to  50  and  the  arrows  represent  the  principal  stresses.  Converging  pairs  represent 
compressible  stresses  and  diverging  pairs  represent  tension  stresses. 


improves  with  increased  number  of  viable  individuals.  In  the  Genetic  Algorithm, 
the  individuals  are  not  viable,  for  instance,  either  because  of  connectivity  prob¬ 
lems,  or  the  removal  of  loaded  elements,  whereas  these  defects  are  prevented  in 
bioTOM.  For  bit-array  representations  that  ratio  can  be  rather  low  as  shown  in  the 
results  for  [14]  whereas  for  the  bioTOM  all  the  individuals  are  feasible.  Another 
performance  measurement  is  the  number  of  FE  evaluations  necessary  until  conver¬ 
gence  is  achieved.  For  this  problem,  the  number  of  FE  evaluations  necessary  for 
the  bioTOM  is  at  lest  two  orders  of  magnitudes  lower  than  for  the  other  methods. 


feasible 

mesh 

FE  evaluations 

bioTOM 

0.260 

100% 

N/A1 

-  300 

Wang  et  al.  [14] 

0.325 

17% 

20  x  10 

-  20000 

Balamurugan  et  al.  [15] 

0.340 

N/A2 

20  x  10 

~  100000 

Hamda  et  al.  [13] 

0.330 

N/A2 

20  x  10 

~  16000 

1Not  Applicable 
2Not  Available 


Table  1:  Comparative  analysis  of  several  techniques  for  the  minimum  weight  prob¬ 
lem. 
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2.2.  MAV  wing 

In  this  second  test  case,  the  problem  of  dual  optimization  of  mass  and  perfor¬ 
mance  of  a  flexible  MAV  wing  is  studied.  The  wing  is  modeled  as  a  membrane  sur¬ 
face  (latex  rubber)  reinforced  with  a  stiffer  material  (carbon  fiber  laminates) — the 
elastic  properties  for  the  wing  materials  are  given  in  table  2.2.  The  latex  membrane 
determines  the  wing  planform  shape  and  the  carbon  reinforces  the  structure.  The 
goal  is  to  simultaneously  determine  the  wing’s  shapes  and  structure  layouts  with 
best  mass  and  aerodynamic  performance;  the  latter  measured  by  the  wing  finesse: 
L/D  =  Cl/Cd,  where  Cl  and  Co  are  the  lift  and  drag  coefficients,  respectively. 


Rubber 

Carbon 

Young  Modulus  [GPa] 

2  x  1(T3 

317.2 

Poisson  Coefficient 

0.5 

0.31 

Table  2:  Properties  for  the  materials  used  on  the  wings 

To  better  assess  bioTOMagainst  existing  design  methods,  the  present  test  case 
is  divided  into  four  parts.  In  the  first  three,  the  shape  is  fixed  and  only  the  topology 
is  optimized.  This  allowed  for  a  direct  comparison  between  bioTOMand  current 
popular  methodologies  for  topology  optimization.  Thus  a  rectangular  planform  is 
chosen  with  a  chord  equal  to  c  =  0.15  cm  and  a  span  equal  to  b  =  0.30  cm— 
corresponding  to  an  area  of  A  =  0.045  m2  and  an  aspect  ratio  (AR  =  ( b)2/A )  of 
2.  In  the  fourth  part,  the  versatility  of  bioTOM  is  explored  by  allowing  wing’s 
planform  shape  to  vary,  while  keeping  the  wing  area  constant.  For  all  cases  the 
undeformed  shape  of  the  wing  is  a  flat  plate  and  the  final  wing  shape  is  solely  the 
result  of  the  fluid-structure  interaction.,  m,  m2 

The  wing  is  positioned  relative  to  the  incoming  flow  at  an  angle  of  attack  (a) 
that  is  assumed  to  be  fixed  and  equal  to  3°  for  all  cases. 

The  symmetry  of  the  problem  allows  to  reduce  the  computational  mesh  by 
simulating  only  one  side  of  the  wing.  At  the  symmetry  axis,  which  corresponds 
to  the  wing’s  root,  the  structure  is  assumed  to  be  clamped,  i.e.,  zero  deformation 
and  zero  angular  deformation.  Figure  5  depicts  the  problem  geometry  for  the 
rectangular  wing  (a)  and  for  the  variable  shape  wing  (b).  The  flow  conditions 
correspond  to  a  Reynolds  number  of  105,  a  typical  value  for  these  kind  of  aerial 
vehicles  [16].  The  Reynolds  is  defined  as  Re  =  where  Lc  is  a  characteristic 

length  for  the  problem.  For  all  cases  Lc  is  set  to  be  equal  to  the  chord  of  the 
rectangular  wing,  Lc  =  c  =  0.15  rn.  The  air  properties  are  given  below  in  the 
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Figure  5:  The  geometry  of  the  problem:  (a)  rectangular  wing  and  (b)  the  wing  of 
variable  shape.  Both  wings  have  the  same  planform  area.  The  flow  conditions  are 
the  same  for  both  cases.  The  computational  domains  are  highlighted  in  red. 

table  2.2  and  correspond  to  the  air  properties  for  the  standard  atmosphere  at  sea 
level.  From  these  values  the  velocity  of  the  flow  can  be  readily  calculated  to  be 
V  =  13  ms-1. 


Viscosity  [kg/sm]  1.785  x  10  5 
Density  [ kgm~ 3]  1.225 


Table  3:  Air  properties  at  sea  level  for  the  standard  atmosphere. 


As  previously  explained,  the  topology  generated  by  the  Map  L-system  must 
be  given  a  physical  interpretation.  In  this  test  case,  the  edges  represent  carbon 
elements  and  the  interior  regions  of  the  cells  represent  the  latex  membrane. 

The  relevant  parameters  for  the  Map  L-system  follow  next. 

•  The  initial  map  is  bounded  by  the  rectangle  [0,  0.15]  x  [0,  0.15]  which  corre¬ 
sponds  to  one  half  of  the  rectangular  wing. 

•  The  number  of  developmental  stages,  N Levels  is  an  optimization  parameter 
in  the  range  [3,  6]. 

•  The  alphabet  E  contains  8  letters  excluding  the  terminal  symbol  x. 

•  The  production  rules  have  8  tokens. 
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•  The  special  symbol  *  and  /  are  used  to  control  the  edges  thickness. 

•  The  equilibrium  computation  is  bypassed. 

For  the  fourth  problem,  the  wing  shape  is  also  a  design  parameter  and,  as 
such,  it  must  be  encoded  in  the  individual’s  genome.  In  this  work  the  shape  of  the 
wing  is  constructed  from  2N  parameters  appended  to  the  genome.  The  first  N 
elements  consist  on  the  coordinates  (x,y)  that  define  the  wing  leading  edge  (L.E.) 
P  =  {Pi,--  -  ,  Pat}  -  see  figure  7  where  N  =  5.  The  other  N  elements  correspond 
to  the  wing  chords  C  =  {ci,  •  •  •  ,  cjv}  at  the  spanwise  coordinates  given  by  the  x 
component  of  P.  The  trailing  edge  (T.E.)  chordwise  coordinates  (y  direction)  are 
simply  given  by 

Qy  =  {Ply  +  Cl,  '  •  •  ,  -PjVy  +  Cjv}  (3) 

and  the  x  coordinates  are  Qx  =  Px.  The  leading  edge  and  trailing  edge  coordinates 
are  then  connected  with  splines  to  ensure  smooth  curves  and  these  two  curves  are 
connected  at  the  root  and  tip  by  straight  lines,  forming  the  closed  curve  S.  This 
method  is  unbiased  toward  any  particular  kind  of  wing  planform,  and  can  generate 
a  large  variability  of  wing  shapes  as  figure  6  demonstrates. 


Figure  6:  Examples  of  several  wing  planforms  generated  randomly. 

With  the  planform  shape  S  determined,  all  the  vertices  in  the  Map  L-system 
topology  are  linearly  interpolated  such  that  the  square  [0,  0.15]  x  [0,  0.15]  is  trans¬ 
formed  into  the  shape  S.  The  outcome  is  a  smooth  and  complex  topology  as  shown 
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in  figure  7(c).  Since  after  this  transformation  the  wing  area  and  JR  are  unspecified 
values  it  is  necessary  to  scale  the  geometry  in  order  to  obtain  the  final  topology. 
The  wing  JR  is  also  an  optimization  parameter  encoded  in  the  individuals  genotype 
and  it  is  always  within  the  range  0.5  <  JR  <  10. 

To  obtain  the  desired  values  for  area  and  JR  the  scaling  factors  in  the  chordwise 
direction,  Xsf,  and  in  the  spanwise  direction,  Ysf  are  in  general  distinct,  and  can 
be  calculated  as: 


Xsf  = 

A 

A0Ysf 

b 

Ysf  = 

bo 

b  = 

Vjra 

where  b0  and  b  are  the  wingspans  for  before  and  after  the  scaling,  respectively. 
A0  and  A  are  the  wing  areas  for  before  and  after  the  scaling,  respectively.  This 
process  is  exemplified  in  figure  7 

With  the  geometry  so  defined,  it  is  possible  now  to  build  the  structure  model 
and  determine  the  wing  deformation  using  the  FE  method  for  a  specific  wing  load. 
As  mentioned  above  the  edges  of  the  Map  L-system  topology  are  interpreted  as 
Euler  beams  with  square  cross  section  made  out  of  carbon  laminate.  The  beams 
are  characterized  by  the  material  properties  Ec,  uc,  pc  and  cross  section  properties: 
area  Ac  and  moments  of  area  Iy,  Iz  and  J .  As  opposed  to  the  cantilever  problem 
shown  before,  here  the  inner  region  of  the  Map  L-system  has  a  physical  meaning: 
the  different  cells  represent  latex  shells  characterized  by  the  material  properties 
Em,  vm,  pm  and  thickness  tm.  The  shell  thickness  is  the  same  for  all  the  four 
optimizations,  tm  =  0.1  mm,  but  the  square  cross  section  of  the  beams  is  not  fixed. 
As  explained  previously,  edge  properties  can  also  be  subjected  to  the  Map  L-system 
production  rules.  That  feature  is  used  in  this  case  to  establish  the  dimensions  of 
the  square  cross  section  of  each  beam.  In  this  case  the  Map  L-system  returns  for 
each  edge  the  property  thicknessRatio  which  is  a  value  in  the  range  [1/3,  3].  That 
value  is  used  to  determined  the  actual  beam  thickness  given  by  t0  x  thicknessRatio, 
where  t0  =  1  mm.  Thus  the  beam  thickness  varies  in  the  range  [0.33,  3]  mm. 

The  geometry  modeling,  meshing,  solving  and  postprocessing  are  done  through 
scripting  in  COMSOL  allowing  for  the  automated  study  of  the  wing  designs.  A 
study  was  performed  to  determine  the  mesh  density  necessary  to  attain  mesh 
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Figure  7:  Example  that  shows  how  the  wing  of  variable  shape  is  ’’build”  from 
the  Map  L-system  topology  and  the  shape  S.( a)  Map  L  systems  topology,  (b) 
wing  platform  shape  S,  (c)  wing  showing  the  internal  structure  and  (d)  the  final 
topology  with  the  correct  area  and  At,  after  applying  the  scaling. 


independence — see  Figure  8.T  ypically,  in  a  2GHz,  4.00  GB  computer  the  FE 
method  solution  is  calculated  after  2  seconds. 

The  procedure  explained  above  is  independent  of  the  wing  load.  To  obtain 
the  steady  state  deformation  for  the  wing  the  aerodynamic  pressure  loads  must 
be  calculated.  There  are  many  computational  tools  that  can  be  used  for  this 
calculation,  such  as,  commercial  or  open  source  softwares  that  solve  the  governing 
equations  -  conservation  of  mass  and  conservation  of  momentum,  that  govern  the 
incompressible  viscous  flows.  These  tools  can  be  very  accurate  but  are  in  general 
computationally  demanding  due  to  the  nonlinear  nature  of  the  conservation  of 
momentum  equations  and  the  necessity  of  dense  meshes  in  the  domain  surrounding 
the  wing.  Since  the  bioTOM  requires  thousands  of  these  computations  these  tools 
can  not  be  used  due  to  the  prohibiting  computational  times.  As  an  inexpensive 
surrogate  for  the  more  accurate  tools  a  Vortex  Lattice  Method  (VLM)  is  used 
here.  The  computational  effort  required  by  the  vortex  lattice  method  is  almost 
negligible  when  compared  to  the  codes  mentioned  above  given  that  it  requires 
only  the  meshing  of  the  wing  surface  and  the  solution  of  a  single  linear  system  of 
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Figure  8:  (a)  An  example  of  a  COMSOL  mesh,  (b)  An  example  of  a  COMSOL 
geometry.  The  colors  represent  the  edge  thickness.  Dark  red:  thick  beams,  dark 
blue:  thin  beams,  (c)  FE  solution.  The  colormap  depicts  the  wing  deformation  in 
m  for  a  constant  load  of  10  Pa 


equations.  So  that  the  velocities  are  given  by  the  gradient  of  the  potential  0: 

v  =  V0,  (5) 

where  the  velocity  potential  is  the  solution  of  the  Laplace  equation 

V20  =  O,  (6) 


and  the  boundary  conditions  are:  velocity  must  tend  to  the  approaching  velocity 
in  the  farfield 

v  — >  {V cos(a),  V sin(«)},  (7) 

and  the  velocity  normal  to  the  wing  surface  must  be  zero. 


dv 

dn 


=  0. 


(8) 
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The  pressure  can  be  computed  from  the  Bernoulli  equation: 

m/2  m/2 


where  V  =  |v 

The  denomination  of  vortex  lattice  method  encompasses  a  myriad  of  methods. 
This  work  uses  a  linear  vortex  lattice  method  as  described  in  [17].  The  vortex 
lattice  method  approximates  the  continuous  distribution  of  bound  vorticity  by 
discretizing  the  wing  into  a  paneled  grid  of  N  spanwise  panels  and  M  chordwise 
panels  as  exemplified  in  figure  9,  and  placing  a  vortex  ring  upon  each  panel.  The 
vortex  ring  is  created  by  four  vortex  filaments  located  at  the  edges  of  the  panel. 
Each  vortex  filament  creates  a  velocity  whose  magnitude  is  assumed  to  be  governed 


Figure  9:  (a)  Example  of  the  paneling  of  a  generic  wing  with  N  =  5  and  M  =  3. 
The  black  dots  are  the  wing  grid,  and  the  red  dots  are  the  points  where  the  vortex 
is  placed  -  collocation  points.  Three  panels  are  show,  (b)  shows  a  typical  panel 
for  which  i  =  2  and  j  —  3.  n,-?  is  the  normal  vector  to  the  panel  ij  and  F^  is  the 
vortex  strength  [m2/s]. 

by  the  Biot-Savart  law  [17].  Furthermore,  a  control  point  is  placed  at  the  three- 
quarter-chord  point  of  each  panel.  The  tangency  condition  is  applied  (i.e.,  the 
wing  becomes  a  streamline  of  the  flow)  by  stipulating  that  the  induced  flow  (from 
the  vortex  rings)  along  the  outward  normal  at  each  control  point  exactly  cancels 
with  that  caused  by  the  free-stream  velocity.  The  following  system  of  equations 
results: 


L  Oil 

0-12  '  ' 

Cilm  ^ 

frA 

/  (V cos(a),  0,  V sin(a))  •  ni  \ 

021 

022 

^2  m 

r2 

—  — 

( V cos(cr),  0,  V  sin(tt))  •  n2 

y(2mi 

0-rn2 

&  mm  J 

\W 

\(y  cos(«),  0,  V sin(cr))  •  n mJ 
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r,(  is  the  unknown  circulation  strength  of  each  vortex  ring,  and  m  =  N  x  M.  In 
order  to  obtain  this  system  of  equations  the  panels  must  be  relabeled  in  a  sequential 
manner  from  0  to  m  =  N  x  M.  For  the  structured  grid  that  is  easily  done  with 
the  aid  of  a  variable  K  =  (j  —  1)  x  N  +  i.  As  an  example,  the  panel  in  figure  9(b) 
is  relabeled  from  (2,  3)  to  12.  The  elements  in  the  matrix  from  equation  10  are  the 
influence  coefficients,  which  are  given  by: 

aij  =  (u0,v0,  w0)ij  ■  ni,  (11) 

where  (u0,v0,w0)ij  are  the  velocities  induced  at  the  i-th  collocation  point  by  the 
j-th.  The  indices  i  and  j  are  sequential  indices,  each  one  identifies  a  different 
panel.  The  solution  of  the  linear  system  (10)  yields  the  strength  circulation  for 
each  panel.  Further  details  and  information  can  be  found  in  the  [17].  Before  using 
the  code  for  the  optimization,  the  code  was  thoroughly  verified  against  analytical 
solutions  and  reference  numerical  data  [17]. 

At  this  point  the  means  to  calculate  the  wing  deformation  for  any  given  out 
of  plane  load  and  the  means  to  determine  the  aerodynamic  load  for  any  given 
wing  geometry  are  defined.  The  following  step  is  to  incorporate  these  two  solvers 
in  a  fluid-structure  interaction  solver  in  order  to  obtain  the  steady  state  wing 
deformation.  Given  that  the  FE  method  here  used  is  the  one  implemented  in 
the  commercial  package  COMSOL,  and  the  vortex  lattice  method  is  coded  and 
compiled  in  FORTRAN,  a  segregated  solver  was  used. 

The  segregated  solver  was  implemented  in  MATLAB  and  consists  of  a  fixed 
point  iteration  with  intermediate  update.  Thus  the  aerodynamic  load  is  computed 
first  for  the  deformed  wing  and  then  the  wing  deformation  is  calculated  for  that 
load  and  so  on.  The  communication  between  the  FE  method  and  the  vortex  lattice 
method  is  processed  through  ASCII  files. 

The  loop  continues  until  the  stopping  criteria  are  satisfied.  In  this  case  the 
solver  stops  if  the  value  of  the  finesse  Cl/Cd  converges  or  if  the  maximum  num¬ 
ber  of  iterations  is  achieved.  The  last  condition  is  necessary  to  handle  divergence 
of  the  segregated  approach,  which  was  initially  observed  for  the  most  compliant 
structures — the  large  deformations  lead  frequently  to  unstable  solutions.  For  most 
of  these  cases  a  stable  solution  can  be  obtained  if  relaxation  is  used:  the  current 
deformation  is  a  linear  interpolation  between  the  deformation  from  the  previous  it¬ 
eration  and  the  one  from  the  current  iteration  u  =  r*uNew  +  (l-r)*u01d.  Where 
r  is  the  relaxation  factor  that  varies  between  0  and  1. 

In  pseudo  code  the  aeroelastic  solver  can  be  written  as: 

ContinueLoop  =  true 
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uNew  =  0; 
uOld  =  0; 
u  =  0; 
it  =  1 

while  ContinueLoop: 

P  =  computeWingLoad(u) 

uNew  =  computeWIngDef ormation(P) 

if  StopLoop  ==  true: 

ContinueLoop  =  false; 
else : 

ContinueLoop  =  true; 
endif 


it=it+l ; 

u  =  r*uNew  +  (l-r)*u01d; 
uOld  =  uNew; 
end  while 

where  the  relevant  subroutines  and  variables  are  listed  next. 

•  u,  uNew  and  uOld  are  the  current  (relaxed)  deformation,  the  latest  deforma¬ 
tion  and  the  deformation  form  the  privious  iteration,  respectively. 

•  P  is  the  current  aerodynamic  load. 

•  ComputeWingLoad:  returns  the  aerodynamic  load  for  a  given  wing  geometry 
and  deformation. 

•  computeWIngDef ormation:  returns  wing  deformation  uNew  for  a  given  wing 
load  P. 

•  StopLoop:  returns  true  if  the  one  of  the  stopping  criteria  is  satisfied.  Returns 
false  otherwise. 


Just  as  in  the  previous  test  case,  it  is  necessary  to  define  the  objective  vector. 
In  this  test  case,  there  are  two  goals:  to  minimize  the  wing  mass,  and  to  maximize 
the  wing  aerodynamic  performance  as  given  by  the  finesse  Cl/C e>.  The  finesse 
for  the  deformed  wing  is  calculated  as  explained  previously.  The  mass  of  the 
wing  is  calculated  by  integrating  the  density  along  the  beam  cross  section  or  panel 
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thickness  as  appropriated.  Upon  determine  these  two  quantities  the  fitness  of  a 
particular  design  is  determined: 

{F(x),  W (x)}  +  P§  (12) 

where  F  is  the  negative  value  for  the  wing’s  finesse  and  W  is  the  wing’s  weight, 
x  is  the  vector  of  reals  that  encodes  the  genome  and  P§  is  the  penalization  for 
maximum  deformation,  given  by: 

Ps  =  5  (emax(<5  max  $ limt  °)  -  1) 

The  inclusion  of  this  term  is  necessary  given  that  for  individuals  with  less  stiff 
elements  the  deformations  can  be  very  large,  which  is  not  desirable  for  several 
reasons.  In  the  first  place  the  validity  of  the  equilibrium  equations  solved  in  the 
FE  method  might  be  compromised  since  they  follow  from  the  linear  elasticity 
framework  which  assumes  small  deformations.  Secondly  the  wing  with  large  de¬ 
formations  could  provoke  massive  flow  separations  which  would  render  the  vortex 
lattice  method  calculations  invalid  as  remarked  before.  As  for  stress  levels,  they 
are  assumed  to  be  small  such  that  there  is  no  risk  of  material  failure  and  therefore 
no  stress  constraint  or  penalization  needs  to  be  considered. 

Methodology  for  the  Genetic  Algorithm  with  bit-array  representation 

For  this  methodology  the  wing  structure  is  given  by  a  quadrangular  grid  of 
N  x  M  elements,  where  N  and  M  are  the  number  of  panels  in  the  chordwise  and 
spanwise  directions,  respectively.  This  is  the  same  grid  used  for  the  vortex  lattice 
method  which  was  previously  explained.  However,  now  each  grid  element  also 
represents  a  panel  or  shell  as  exemplified  in  figure  11.  To  each  panel  is  assigned  a 
the  value  0  or  1.  The  former  indicates  that  the  material  of  the  panel  is  latex,  and 
the  latter  indicates  the  panel  is  made  of  carbon  laminate. 

The  wing  deformation  is  computed  with  the  FE  method  using  shell  elements, 
whose  material  can  be  latex  or  carbon  laminate.  The  material  properties  are  given 
in  table  2.2.  An  important  difference  between  bioTOM  and  Genetic  Algorithm  is 
that  now  the  carbon  laminate  is  in  the  form  of  panels  and  no  longer  in  the  form  of 
beams.  The  thickness  of  the  carbon  laminate  panel  is  the  same  as  the  latex  panels 
tm  =  0.1  mm. 

The  vortex  lattice  method  is  exactly  as  for  the  bioTOM  as  well  as  the  aeroelastic 
solver. 

The  optimization  is  performed  with  the  Genetic  Algorithm  for  witch  the  genome 
is  a  vector  of  iV  x  M  binary  elements  that  determine  the  wing  structure  as  ex¬ 
plained  above.  Figure  10  exemplifies  the  decoding  for  the  bit-array  representation. 
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Figure  10:  Decoding  step  for  the  bit-array  representation,  (a)  The  genome;  (b)  the 
resulting  wing  structure:  white  squares  represent  latex  panel  and  the  black  squares 
represent  carbon  panels. 

Methodology  for  the  gradient  method  with  SIMP 

Here  the  solid  isotropic  material  with  penalization  (SIMP)  method  is  briefly 
outline.  The  interested  reader  may  consult  [18],  [16].  For  this  method,  the  topology 
of  the  wing  is  defined  in  a  similar  way  to  the  bit-array  representation.  The  main 
difference  is  that  to  each  panel  is  assigned  a  value  that  varies  between  0  and  1. 
This  value,  for  a  given  panel  i,  1  <  i  <  M  x  N,  represents  the  density  variable 
Xi,  where  the  word  density  stands  for  the  relative  percentage  of  each  material. 
For  example,  a  density  value  of  Xi  =  0  means  that  the  panel  i  is  made  out  of 
latex  membrane,  a  density  value  of  X*  =  1  represents  a  panel  made  out  of  carbon 
laminate  and  an  intermediate  value  represents  a  theoretical  inexistent  material 
that  is  a  mix  of  both.  The  efficiency  of  this  method  increases  when  X  is  penalized 
[18].  Penalization  is  accomplished  by  replacing  the  changing  the  density  from  X,t 
to  Xf ,  where  p  >  0  is  the  penalization  power. 

The  usage  of  an  inexistent  material  arises  from  the  necessity  of  having  a  smooth 
function  that  can  be  optimized  using  the  gradient  method.  This  is  a  theoretical 
artifice  that  has  been  widely  and  successfully  used  in  topology  optimization  [18]. 
Formally  the  problem  can  be  stated  as: 

minimize  F(x),x  (13) 

such  that  IT(x)  =  —  Wum 

where  F  and  W  are  the  negative  of  the  finesse  and  the  wing’s  weight  as  defined 
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previously.  In  this  case  the  exact  expressions  for  the  mass  constraint  can  be  easily 
obtained: 

NxM 

m  =  E  mi  (14) 

2=1 

where  nrii  is  the  mass  of  the  i-th  panel  which  is  given  by 

TTii  X{  •  Tflc  +  (1  X^  •  (15) 


where  mc  and  mm  correspond  to  the  mass  of  the  carbon  laminate  panel  and  to  the 
mass  of  the  latex  panel,  respectively.  Replacing  15  in  to  14  yields  a  linear  equation 
for  the  constraint. 

NxM 


m  = 


E  Xi  '  +  (!  “  Xi)mr, 


2=1 


which  can  be  recast  as: 


(11  •••  11) 


X1  \ 
X2 

\XNxM  J 


m  —  N  x  M  • 


(16) 


(17) 


The  optimization  was  done 
finesse. 


keeping  the  mass  constant  and  maximizing  the 


dF 

dX 

For  a  given  X  the  finesses  is  calculated  with  the  areoelastic  solver  explained  previ¬ 
ously.  The  only  difference  is  that  now,  for  the  FE  method  calculation  the  material 
properties  for  the  shell  elements  are  those  of  the  fictitious  material  determined  by 
the  X,t .  In  practical  terms  the  stiffness  matrix  Ke  of  each  finite  element  is  given 
by: 

Ke  =  (Kc  -  Km)  •  XI  +  Km  (19) 

where  Xe  is  the  density  of  the  element,  and  p  is  the  nonlinear  penalization  power. 

With  the  structure  determined  as  exemplified  in  figure  11,  the  wing  deforma¬ 
tion  is  computed  using  the  FE  method  implemented  in  MATLAB1.  The  grid  for  the 
FE  method  consists  in  a  triangular  mesh  that  is  obtained  by  bisecting  each  panel 
diagonally  as  depicted  in  figure  11(c).  For  the  results  displayed  next  N  =  M  =  30, 

1The  Pis  thank  Bret  Stanford  for  providing  the  initial  SIMP  code. 
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which  yields  900  panels.  These  figures  were  dictated  mainly  because  of  compu¬ 
tational  limitations,  ideally  more  elements  would  be  desired  for  reasons  of  grid 
resolution,  since  the  distribution  of  membrane  and  carbon  elements  is  dependent 
on  the  size  of  the  panels. 

Here  it  should  be  pointed  out  a  major  difference  between  the  bioTOM  and  the 
other  two  methodologies:  the  mesh  independence  that  characterizes  the  bioTOM. 
Indeed,  nowhere  in  the  bioTOM  formulation  a  mesh  or  ground  structure  is  nec¬ 
essary.  That  has  a  tremendous  effect  on  the  maximum  resolution — in  the  sense 
of  the  smallest  element  on  the  wing  structure — that  the  structure  can  display. 
For  the  non-bioTOM  approaches  such  as  these  two  the  resolution  is  determined 
by  the  grid,  in  this  case  with  the  mesh  of  30  elements  per  side,  the  resolution  is 
150/30  =  5  mm,  whereas  for  the  bioTOM  that  value  is  only  dependent  on  the 
limits  imposed  by  the  user.  As  explained  above  the  bioTOM  smallest  element  can 
be  0.33  mm.  To  obtain  that  resolution  with  a  mesh  such  as  the  ones  here  used, 
455  elements  per  side  would  be  necessary,  yielding  a  total  of  207,  025  panels.  A 
grid  this  fine  would  be  computational  prohibitive. 


(a)  (b)  (c) 

Figure  11:  Wing  structure  for  the  Optimization  1  (a)  and  Optimization  2  (b). 
The  grey  colormap  represent  the  density  variable  X.  Black  is  for  1  and  white  is  0. 
(c)  depicts  the  triangular  mesh  for  the  FE  method.  In  this  example  N  =  M  =  8. 

In  the  following,  the  results  for  the  four  optimizations  are  presented.  First 
it  is  shown  the  results  obtained  with  the  SIMP  then  the  results  for  the  Genetic 
Algorithm  with  a  bit-array  representation  and  finally,  it  is  presented  the  result  for 
the  bioTOM.  For  reference,  at  this  flight  condition  the  finesse  for  a  rigid  rectangular 
wing  is  equal  to  2.49. 

The  optimization  with  SIMP  ran  for  100  iterations  or  until  the  finesse  con¬ 
verges.  A  typical  convergence  history  is  shown  in  figure  12. 
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Figure  12:  Convergence  history  for  the  finesse  optimization. 


Several  values  for  the  mass  were  studied  spanning  from  a  wing  made  out  of 
latex  to  a  wing  made  out  of  carbon.  For  easy  of  reading,  in  the  results  below  the 
mass  is  represented  by  the  percentage  carbon  laminate  coverage. 

As  explained  previously,  the  output  of  the  minimization  is  the  vector  of  densi¬ 
ties  X  which  requires  post  processing  given  that  values  between  0  and  1  represent 
a  fictitious  material.  This  is  a  well  know  limitation  of  homogeneous  methods  and 
there  are  several  techniques  to  addresses  this  issue  [16].  Here  a  post-processing  is 
employed  by  applying  the  transformation  to  the  elements  Xf 

0  if  Xi  <  e 
1  if  Xi  >  e 

where  e  is  a  parameter  controlled  by  the  user.  In  the  results  that  follow  e  = 
0.5  which  corresponds  to  rounding  the  elements  of  X.  This  correction  alters  the 
structure  and  the  finesse  as  shown  in  figure  13,  but  maintain  the  correct  trend.  In 
this  figure  the  gray  line  refers  to  the  results  without  postprocessing  and  is  shown 
just  for  reference  since  it  represents  fictitious  structures.  The  “true”  results  are 
shown  in  black.  Figure  14  shows  the  best  wing  structure  of  all  cases  studied. 
The  optimum  level  of  carbon  laminate  coverage  was  found  to  be  approximate  30% 
without  postprocessing  -  figure  14(a).  After  the  postprocessing  -  figure  14(a),  the 
percentage  of  carbon  laminate  decreases  to  20%  but  the  finesse  remains  almost  the 
same  at  5.4  as  shown  in  figure  fg:bestSimp. 

This  figure  clearly  shows  that  there  is  a  optimum  level  of  carbon  coverage  such 
that  the  wing  is  neither  to  stiff  nor  to  compliant. 
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Figure  13:  The  finesse  variation  with  the  carbon  laminate  coverage. 
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Flow 


(a)  (b) 

Figure  14:  Best  wing  structure.  Only  half  of  the  wing  is  shown.  The  arrow  shows 
the  flow  direction,  (a)  Without  postprocessing,  (b)  with  postprocessing 

In  this  segment,  the  results  for  the  square  wing  optimized  with  the  Genetic 
Algorithm  and  the  bit-array  representation  for  the  wing’s  structure  are  presented. 
The  Genetic  Algorithm  ran  for  100  generations  of  200  individuals.  Figure  16 
shows  all  the  individuals  calculated  during  the  simulations.  The  Pareto  front  is 
highlighted  in  blue.  The  horizontal  axis  represents  the  designs  finesse  and  the 
vertical  axis  represents  the  wing’s  mass.  Given  that  the  carbon  coverage  is  a 
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(a)  (b) 


Figure  15:  The  wing  deformation  for  the  best  structure,  (a)  The  color  map  shows 
the  wing  deformation  in  m.  The  contours  identify  carbon  laminate  regions. (b) 
Three  different  perspective  views  of  the  wing  deformation.  The  arrows  show  the 
incoming  flow  direction  for  orientation. 
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Figure  16:  Pareto  front  for  Genetic  Algorithm  with  the  bit-array  representation. 

relevant  parameter  as  demonstrated  previously,  that  value  can  be  read  in  the  left 
vertical  axis. 

There  are  nine  designs  in  the  Pareto  front.  Depending  on  design  goals  any  of 
these  individuals  can  be  selected.  Three  of  these  designs  order  form  highest  to 
lowest  finesse  can  be  seen  in  figures  17,  18  and  19. 
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The  results  for  the  multiobjective  bioTOM  optimization  for  the  square  wing 
follow  next.  These  results  were  obtained  with  100  generations  of  200  individuals. 
All  individual  whose  fitness  are  in  the  range  [1,  6]  x  [0.001, 0.007]  are  shown  in  figure 
20.  The  Pareto  front  is  displayed  in  blue.  The  imediate  conclusion  by  comparing 
to  the  results  from  the  two  previous  results  is  that  the  bioTOM  finds  designs  whose 
finesse  is  as  good  as  the  ones  for  the  SIMP,  but  as  the  mass  decreases  bioTOM 
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Figure  20:  Pareto  front  for  Genetic  Algorithm  with  the  bit-array  representation. 

outperforms  SIMP.  Three  different  designs  are  shown  in  figures  21  to  23.  Each 
panel  shows  from  left  to  right:  the  wing  design  with  the  carbon  beams  colored  by 
their  thickness  in  mm,  the  steady-state  wing  deformation  in  m  and  the  pressure 
coefficient.  This  designs  are  striking  for  their  simplicity.  Indeed  the  evolutionary 


0  0.025  0.05 

(b) 

Figure  21:  Pareto  front  design:  finess=5.40,  mass=0.0026. 

algorithm  shows  that  very  few  to  none  beams  are  necessary  within  the  wing.  All 
emphasis  is  in  the  beams  at  the  edge  of  the  wing.  Figure  24  displays  the  structure 
of  the  wing  with  carbon  beams  of  different  thicknesses  for  the  design  with  the  best 
finesse. 
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Figure  22:  Pareto  front  design:  finess=5.26,  mass=0.0025 
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(a)  (b)  (c) 

Figure  23:  Pareto  front  design:  finess=3,53,  mass=0.0023 

Finally  the  results  for  the  wing  with  a  variable  shape  are  shown.  The  optimiza¬ 
tion  was  ran  with  a  population  of  two  hundred  individuals  and  for  one  hundred 
generations.  Figure  25  shows  all  designs  studied  whose  fitness  vector  falls  in  the 
range  [0,10]  x  [0,0.01].  Some  Pareto  designs  are  shown  in  figures  26-29.  Each 
of  these  figures  depicts:  the  (a)  the  wing  geometry  or  skeleton  with  the  beams 
thickness  given  by  the  colormap;  (b)  the  wing  deformation  in  meters ;  and  (c) 
the  pressure  coefficient  distribution.  All  the  figures  assume  the  same  color  maps. 
These  figures  are  shown  in  descending  order  of  finesse. 

The  amplification  factor  is  not  the  same  for  every  figure,  or  otherwise  some 
figures  would  be  very  small  because  the  wing  M  differs  from  design  to  design, 
hence  the  wing  shapes  can  not  be  compared  using  these  figures.  Figure  30  depicts 
the  wing  designs  in  the  same  scale.  The  labels  in  this  figure  show  the  wing  aspect 
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Figure  24:  The  actual  wing  structure  showing  the  carbon  beams  with  different 
thicknesses. 
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Figure  25:  Pareto  front  and  accumulated  population  at  the  one  hundredth  gener¬ 
ation. 

ratio,  M. 

The  layouts  on  the  Pareto  front  conveys  a  trend  for  increased  “area  coverage” 
by  the  reinforcements  as  more  mass  becomes  available.  Thus,  small  mass  Pareto 
designs  tend  to  concentrate  their  limited  mass  along  critical  axes,  whereas  with 
more  mass  available  the  the  number  and  spread  of  reinforcements  increase.  This 
increase  in  area  coverage  with  more  mass  available  can  be  understood  by  noticing 
that,  at  low  mass,  spreading  this  low  amount  of  reinforcement  throughout  the  wing 
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(a)  (b)  (c) 

Figure  26:  Pareto  design:  finess=8.87,  mass=0.0054.  (a)  geometry.  The  color  bar 
represents  the  beams  thickness  in  mm.  (b)  wing  displacement  in  m.  (c)  pressure 
coefficient. 

box  area  would  not  strengthen  any  specific  wing  region.  And  would,  as  a  conse¬ 
quence,  provide  little  stress  relief  on  the  wing  box.  By  contrast,  with  increased 
mass  available,  the  topology  can  strengthen  the  critical  axes  of  the  structure  and 
also  reduce  the  stress  with  reinforcements  that  work  together  to  distribute  the  load 
amon 

3.  Fighter  wing  box 

This  next  test  case  consider  the  problem  of  the  design  of  the  wing  box  of  a 
generic  aircraft  fighter.  As  before,  the  topology  generated  by  the  map  L  systems 
is  a  planar  topology  with  no  canonical  structural  meaning  attached  to  it — in  the 
sequel,  this  topology  is  called  the  structure  skeleton.  To  build  a  structure  model 
from  this  topology,  first  it  is  selected  the  following  basic  structural  elements: 

•  Shear  panels; 

•  Posts  and  caps  for  each  shear  panel. 

•  Upper  and  lower  skin  shells  connected  to  the  shear  panels. 

The  shear  panels  are  structural  elements  that  can  support  shear  stress  and 
extensional  force  along  its  edges.  These  elements  have  negligible  bending  stiffness 
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Figure  27:  Pareto  design:  finess=7.67,  mass=0.0035.  (a)  geometry.  The  color  bar 
represents  the  beam’s  thickness  in  mm.  (b)  wing  displacement  in  m.  (c)  pressure 
coefficient  on  the  wing. 


and  axial  membrane  stiffness.  The  posts  and  caps  are  rod  elements  that  connect 
the  shear  panels  and  shells,  and  they  have  negligible  bending  stiffness.  The  upper 
and  lower  skins  are  modeled  as  isotropic  shells. 

These  basic  structural  elements  are  then  combined  into  a  structural  model  for 
the  wing  box.  The  resulting  wing  box  is  simulated  using  the  finite-element  method 
[19,  20,  21].  The  final  finite  element  model  is  explained  and  illustrated  below  for 
the  wing  box  structure  of  a  generic  aircraft  fighter.  In  the  present  test  case,  because 
of  its  widespread  use  in  the  aeronautical  industry,  the  FE  model  was  created  and 
simulated  using  MSC/NASTRAN™.  This  choice  also  highlights  the  versatility  of 
bioTOM  in  working  with  existing  software,  with  minimum  modification. 

The  specific  steps  leading  to  the  finite  element  model  of  the  wing  box  are 
enumerated  below. 

1.  Mesh  generation.  Given  a  layout,  in  this  first  step,  a  two-dimensional 
mesh  is  generated  that  conforms  with  the  planar  topology — see  figure  31. 
For  that  purpose,  an  advancing  front,  Delaunay  triangulation  method  [19]  is 
used.  The  result  is  a  set  of  two-dimensional  grid  points,  and  the  adjacency 
matrices  for  the  nodes,  edges  and  faces  of  the  discretized  domain.  In  the 
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Figure  28:  Pareto  design:  finess=5.29,  mass=0.0028.  (a)  geometry.  The  color  bar 
represents  the  beam’s  thickness  in  mm.  (b)  wing  displacement  in  m.  (c)  pressure 
coefficient  on  the  wing. 


NASTRAN™  card  deck,  the  “GRID”  points  are  written  in  the  three  dimen¬ 
sional  space  by  symmetrically  placing  the  grid  points  above  and  below  the 
wing  symmetry  plane — the  ^-coordinate  of  each  point  is  obtained  accord¬ 
ing  to  an  interpolation  function  that  accounts  for  the  variation  of  the  wing 
thickness. 

2.  Shear  panels.  After  the  grid  points  have  been  created,  as  a  second  step, 
the  cards  for  the  shear  panels  as  “CSHEAR”  elements  are  written.  These 
elements  are  specified  as  quadrilaterals  using  the  grid  points  generated  in 
the  first  step.  The  thickness  of  each  shear  panel  is  a  function  of  the  edge  in 
the  skeleton  onto  which  the  panel  projects.  In  figure  32  it  is  shown  the  shear 
panels  for  the  mesh  generated  in  the  first  step. 

3.  Posts  and  caps.  In  this  third  step,  the  shear  panels  defined  previously 
are  strengthened  with  posts  and  caps.  The  latter  are  rods  that  are  modeled 
as  “CROD”  finite  elements  in  NASTRAN™.  Each  post  is  placed  at  the 
vertices  of  the  skeleton  and  the  caps  run  at  the  lower  and  upper  edges  of 
the  shear  panels.  The  radius  of  each  cap  is  the  same  for  each  edge  of  the 
skeleton.  The  elements  thus  generated  are  illustrated  in  figure  33. 
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Figure  29:  Pareto  design:  finess=3.48,  mass=0.0021.  (a)  geometry.  The  color  bar 
represents  the  beam’s  thickness  in  mm.  (b)  wing  displacement  in  m.  (c)  pressure 
coefficient  on  the  wing. 


Figure  30:  Individuals  from  the  Pareto  front.  The  label  is  the  Ai  of  the  wing. 


4.  Upper  and  lower  skins.  Finally,  the  wing  box  is  closed  with  upper 
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Figure  31:  Skeleton  (left)  of  wing  box  design  and  the  generated  mesh  (right). 


and  lower  skins.  Both  skins  are  modeled  by  triangular  shell  elements  of 
“CTRIA3”  type  in  NASTRAN™  with  equal  in-  and  out-plane  stiffness.  The 
complete  structural  model  is  depicted  in  figure  34. 

With  the  structural  model,  the  analysis  can  be  performed  using  NASTRAN™ 
to  obtain  the  desired  information  on  deformation,  stress  and  mass  of  the  structure. 
Figure  35  shows  the  simulation  result  for  the  static  aeroelastic  problem  of  the  trim 
problem  in  a  six-gravity  maneuver — see  §??  for  details  on  the  physical  problem. 

The  geometry  of  the  problem  consist  of  a  sweep  and  tapered  wing  with  an 
internal  polygonal  area  whose  boundary  the  structural  elements  cannot  cross — 
see  figure  36.  In  this  problem,  the  polygon  area  inside  the  planform  represents 
the  area  that  delimits  a  generic  aircraft  subsystem.  All  components  of  the  wing 
are  build  with  the  same  material  with  assumed  properties:  E  =  lx  107  psi  and 
v  =  0.33.  The  flight  condition  corresponds  to  the  trim  condition  for  a  vertical,  six- 
gravity  maneuver  at  sea-level  and  Mach  number  M  =  0.7.  The  dynamic  pressure 
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Figure  32:  Shear  panel  elements  in  the  wing  box  model. 


is  q  =  20  psi.  The  aerodynamic  load  is  calculated  using  the  Doublet-Latice  Method 
for  interfering  lifting  surfaces  in  subsonic  flow  [22,  23].  Besides  the  aerodynamic 
loads,  a  concentrated  mass  of  2000  I bm  at  x  =  94.15  in,  y  =  47.05  in,  z  =  —2.56  in, 
a  concentrated  mass  at  the  support  point  at  the  root:  x  =  40  in ,  y  =  0  in ,  ^  =  Oin 
and  a  concentrated  load  of  1000  Ibf  at  the  positive  Oz-direction  at  the  tip  point 
x  =  135.31  in,  y  =  138.20  in,  z  =  1.51  in  are  considered.  The  distributed  body 
gravitational  force  field  is  also  accounted  for  and  an  inertial,  non-structural  added 
mass  distributed  over  the  wing  is  used  to  account  for  inertia  effects:  the  “inertia 
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Figure  33:  Post  and  cap  elements  in  the  wing  box  model. 


relief”  from  NATSRAN™. 

The  optimization  was  run  with  a  population  of  two  hundred  individuals  and 
for  fifty  generations.  In  the  runs,  rules  defined  by  six  token  slots,  and  symbols 
controlling  each  of  the  following  structure  properties:  post  radius,  cap  radius, 
skin  thickness  and  shear  panel  thickness  are  used.  The  vector  of  design  variables 
comprised  the  mass  of  the  structure  and  the  maximum  stress  on  the  structure. 
The  latter  is  calculated  among  the  maximum  von  Mises  stress  in  the  skins,  the 
maximum  axial  stress  on  the  rods  and  the  maximum  shear  stress  in  the  shear 
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Figure  34:  Upper  and  lower  skins  and  the  final  wing  box  model. 


panels. 

Figure  37  shows  the  Pareto  front  and  all  individuals  calculated  in  the  fifty 
generations,  whose  normalized  fitness  vector  fall  in  the  unit  [0,  l]2  square.  The 
normalization  is  effected  using  a  reference  mass  of  1500 1  bm  and  the  allowable 
shear  stress  for  aluminum  12,  000  Ibm. 

A  sample  of  Pareto  designs  is  shown  in  figures  38-45.  Each  of  these  figures 
depicts  the  skeleton  on  the  left  and  the  von  Mises  stresses  on  the  deformed  structure 
at  the  right.  The  displacement  is  amplified  for  easy  of  reading — note  that  each 


Figure  35:  Von  Mises  stress  distribution  for  the  wing  box  model. 


figure  assumes  a  different  color  map  and  different  amplification  ratio.  These  figures 
are  shown  in  descending  order  of  mass,  and  both  the  normalized  mass  and  the 
normalized  stress  are  shown  in  the  caption  as  a  percentage  of  their  respective 
normalizing  factors. 

The  layouts  on  the  Pareto  front  conveys  again  a  trend  for  increased  “area  cov¬ 
erage”  by  the  reinforcements  as  more  mass  becomes  available.  Thus,  small  mass 
Pareto  designs  tend  to  concentrate  their  limited  mass  along  critical  axes,  whereas 
with  more  mass  available  the  number  and  spread  of  reinforcements  increase.  This 
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Figure  36:  Schematic  of  the  wing-box  horizontal  projection.  Coordinates 
(x,y)  of  the  geometry  vertices:  1:  (25.0000,0.00000),  2:  (123.222,0.00000), 

3:  (135.313,138.200),  4:  (121.260,138.200),  5:  (93.2080,15.7500), 

6:  (101.877,25.7500),  7:  (106.570,59.3750),  8:  (96.0300,65.0000),  9: 

(85.4890,59.3750),  10:  (80.7960,25.7500).  Length  in  inches. 


increase  in  area  coverage  with  more  mass  available  can  again  be  understood  by 
noticing  that,  at  low  mass,  spreading  this  low  amount  of  reinforcement  through¬ 
out  the  wing  box  area  would  not  strengthen  any  specific  wing  region.  And  would, 
as  a  consequence,  provide  little  stress  relief  on  the  wing  box.  By  contrast,  with 
increased  mass  available,  the  topology  can  strengthen  the  critical  axes  of  the  struc¬ 
ture  and  also  reduce  the  stress  with  reinforcements  that  work  together  to  distribute 
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Figure  37:  Pareto  front  and  accumulated  population  at  the  fiftieth  generation. 


Figure  38:  Pareto  design:  mass  =  95.4%,  stress  =  17.2%. 
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Figure  39:  Pareto  design:  mass  =  41.1%,  stress  =  20.5%. 


Figure  40:  Pareto  design:  mass  =  36.6%,  stress  =  27.1%. 


the  load  among  several  members. 
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Figure  41:  Pareto  design:  mass  =  35.6%,  stress  =  28.5%. 


Figure  42:  Pareto  design:  mass  =  34.9%,  stress  =  35.2%. 


This  increased  area  coverage  demonstrates,  in  passing,  the  gains  of  topologi¬ 
cal  optimization  when  compared  to  sizing  optimization.  Indeed,  with  more  mass 
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Figure  43:  Pareto  design:  mass  =  22.7%,  stress  =  43.0%. 


Figure  44:  Pareto  design:  mass  =  21.8%,  stress  =  58.9%. 


available,  it  is  preferable  to  redistribute  the  reinforcements  among  several  shear 
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Figure  45:  Pareto  design:  mass  =  16.5%,  stress  =  82.9%. 

panels  rather  than  increase  the  thickness  of  the  existing  layouts  of  the  low  mass 
designs. 

The  figures  also  reveal  a  larger  accumulation  of  structural  reinforcements  at 
the  leading-edge  section  of  the  wing,  specially  at  the  root  section.  This  structural 
buildup  at  the  leading-edge/root  is  essential  to  strengthen  the  wing  against  the 
aerodynamic  pressure  loads,  which  peak  at  the  stagnation  fine  close  to  the  leading 
edge. 

4.  Satellite  adapter 

The  last  test  case  in  this  summary,  is  the  design  of  a  structure  for  a  Payload 
Adapter  and  Deployer,  PAD  for  short.  The  PAD’s  function  is  to  accommodate 
multiple  payloads  into  the  payload  bay.  In  this  case  the  payloads  consist  in  one 
large  satellite,  and  six  small  satellites.  The  PAD  with  the  satellites  assembled  is 
shown  in  figure  46. 

Because  the  highest  loads  the  PAD  will  experience  occur  during  launching  only 
that  stage  is  considered — after  the  payloads  detached  from  the  PAD  its  mission 
has  been  accomplished  and  it  is  discarded.  The  exact  loads  experienced  by  the 
structure  are  unknown,  thus  the  conditions  for  this  optimization  correspond  to 
a  worst  case  scenario  in  which  the  assembled  structure — PAD  plus  the  payloads, 
experiences  vertical  and  lateral  accelerations  13  times  grater  than  g. 
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Figure  46:  The  assembled  structure.  1:  Large  satellite;  2:  small  satellites;  3:  the 
PAD  showing  an  initial  design  not  optimized. 


Figure  46  shows  an  initial  design2  for  the  PAD.  This  design  was  determined 
using  just  the  designers  experience,  and  no  automatic  optimization  was  performed, 
apart  from  sizing  by  trial  and  error.  This  design  weights  22  kg  and  satisfies  all  the 
constraints  and  requirements. 

The  goal  for  present  optimization  is  to  reduce  the  mass  of  the  PAD  with¬ 
out  weakening  the  structure  so  that  it  can  withstand  the  launching  phase.  For 
aerospace  applications,  weight  reduction  is  extremely  important  for  several  rea¬ 
sons.  In  first  place,  cost  reduction  as  at  current  prices  hauling  one  kilogram  into 
a  low  earth  orbit  cost  several  thousand  dollars.  Secondly,  the  reduction  of  any 
unnecessary  weight  allows  for  the  inclusion  of  other  payloads  that  can  enhance  the 
mission,  such  as  batteries,  scientific  instruments,  etc. 


2The  Pis  are  grateful  to  Carol  Hude  for  providing  the  initial  design. 
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The  PAD  is  built  from  aerospace  aluminum  with  the  following  material  prop¬ 
erties: 


Density  2800  kg/m 3 

Young’s  Modulus  71  GPa 
Yield  stress  572  MPa 


Table  4:  Material  properties 

There  are  several  elements  of  the  design  are  determined  a  priori  .  For  example, 
due  to  payload  bay  space  constraints  the  geometric  envelop  for  the  PAD  is  well 
defined — see  figure  47.  Also  the  positions  of  certain  elements  are  constrained 
such  that  the  satellites  can  be  attached  to  the  PAD.  These  elements  are  shown  in 
figure  47  and  are  designated  hereafter  as  the  PAD’s  “bare  bones” . 


Figure  47:  The  bare  bones  of  the  PAD.  The  elements  onto  which  the  satellites  will 
be  attached  are  highlighted  in  color,  red:  large  satellite,  blue:  six  small  satellites. 


The  large  void  areas  in  the  three  platforms  (highlighted  in  yellow  in  figure  47) 
are  free  to  be  “filled”  with  topologies  produced  by  the  Map  L- system.  The  opti¬ 
mization  searches  for  the  best  structure  generated  by  the  PAD’s  bare  bones  and 
Map  L-system.  Furthermore,  the  thickness  of  every  beam  in  the  structure  is  also 
a  optimization  variable. 

Other  than  these  geometrical  constraints  that  can  be  easily  fulfilled  during 
the  geometry  modeling,  there  are  maximum  stress  and  maximum  displacement 
constraints.  The  former  must  be  satisfied  to  guarantee  that  the  structure  does  not 
fail — the  latter  assures  that  the  deformed  PAD  will  not  interfere  with  the  small 
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satellites  and  possibly  damage  them,  given  that  upon  assembling  gaps  between 
PAD  and  satellite  are  narrow.  Formally,  this  problem  can  be  stated  as: 


minimize  VF(x)  (21) 

such  that  amax{x)  -  aHm  <  0 
and  8max(pc)  8 nm  A  0 

where  IF(x)  is  the  PAD’s  weight,  crmax(x)  is  the  maximum  stress  experienced 
by  the  structure,  aum  is  maximum  allowed  stress  level,  8max(x)  is  the  maximum 
displacement  in  the  structure  and  8um  is  the  maximum  allowed  displacement.  For 
this  problem  the  maximum  allowed  values  are: 


® lim 


Jlim 


&yield 

U 

2.0  cm 


(22) 


where  n  is  a  safety  factor  equal  to  1.5,  and  the  can  be  found  in  table  4. 

The  PAD’s  structure  can  be  seen  as  the  assembly  of  many  beams  of  variable 
thickness.  These  structures  can  be  modeled  accurately  with  tridimensional  finite 
elements,  however  mesh  generation  can  be  complex  and  may  require  supervision 
precluding  the  sought  after  automatization  of  the  calculus  of  the  fitness  function. 
Thus,  as  in  the  previous  test  cases,  a  simpler  and  less  computational  expensive 
surrogate  model  is  used  to  determine  the  fitness  value.  For  the  surrogate  model 
all  the  tridimensional  beams  are  replaced  by  Euler  beams  that  assume  a  linear 
variation  of  stress  in  the  transversal  direction  with  the  maximum  stress  located  at 
the  external  faces  of  the  beam.  This  approximation  is  valid  because  all  elements  in 
the  PAD  have  a  high  ratio  of  length  to  thickness,  and  only  small  deformations  are 
considered.  The  usage  of  Euler  beams  simplifies  the  mesh  generation  greatly  since 
all  geometric  elements  become  lines,  and  reduces  the  solution  time  enormously. 
The  differences  in  the  thicknesses  of  the  beams  are  taken  in  to  account  trough  the 
area  properties  -  area  and  area  moments  of  inertia. 

5.  Map  L-system  topology 

The  PAD  structure  is  formed  by  the  PAD’s  bare  bones  plus  five  different  topolo¬ 
gies  determined  by  the  Map  L-system  which  occupy  the  voided  areas  in  the  three 
platforms.  The  way  the  five  maps  -  one  for  the  top  platform,  and  two  for  each  of 
the  other  two,  are  assembled  onto  the  PAD’s  bare  bones  is  shown  in  figure  48. 

The  different  Map  L-system  settings  used  in  the  generation  of  these  topologies 
are  listed  bellow. 
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Figure  48:  This  figure  shows  how  the  five  different  map  L-systems  are  assembled 
in  the  PAD  bare  bones 


•  Alphabet  S  has  6  letters  excluding  x. 

•  The  production  rules  have  6  tokens. 

•  The  number  of  developmental  stages  varies  between  3  and  6. 

•  The  equilibrium  calculation  is  bypassed. 

•  Two  of  the  initial  maps  have  6  edges  and  the  other  three  have  4  edges.  The 
coordinates  of  the  vertices  in  the  initial  maps  are  such  that  the  Map  L-system 
topology  fits  in  the  PAD’s  bare  bones  as  shown  in  figure  48. 
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The  PAD  must  withstand  a  vertical  and  lateral  acceleration  13  times  greater 
than  g  with  the  satellites  attached  onto  it.  For  simplicity  each  the  satellite  is 
replaced  by  a  lateral  and  a  vertical  force  of  magnitude  13 g.m  where  rn  is  the 
mass  of  the  satellite.  The  mass  of  the  large  satellite  is  mi  =  50  kg  and  the  mass 
of  each  one  of  the  small  satellites  is  m2  =  5.25  kg.  These  forces  are  applied  to 
correspondent  satellite  center  of  gravity  as  shown  in  figure  49.  These  forces  are 
then  transfered  to  the  PAD  as  distributed  forces  (force  per  unit  of  length)  and 
moments.  Furthermore  body  forces  due  to  the  structure’s  own  weight  are  also 
considered. 


Figure  49:  Loads  due  to  the  satellites  applied  at  the  respective  centers  of  gravity. 


The  calculation  of  this  forces  and  moments  due  to  the  satellites  weight  is  done 
in  a  simple  way.  For  example,  all  the  effects  caused  by  the  presence  of  the  large 
satellite  will  be  transfered  to  the  its  attachment  points,  in  this  case  the  ring  in 
the  top  platform.  The  vertical  acceleration  experienced  by  the  satellite  causes  a 
distributed  load  in  the  vertical  direction,  whose  magnitude  is  calculated  by  dividing 
the  satellite’s  weight  at  13g’s  by  the  perimeter  of  the  ring.  As  for  the  lateral 
acceleration  experienced  by  the  satellite  its  effect  is  twofold.  In  first  place,  it  causes 
a  lateral  distributed  force  whose  direction  is  the  same  as  the  lateral  acceleration 
and  whose  magnitude  is  same  as  for  the  vertical  force.  Secondly,  it  also  provokes  a 
twisting  moment  on  the  ring.  This  moment  is  simulated  by  applying  a  distributed 
force  onto  the  ring  as  shown  in  figure  50.  The  magnitude  of  this  force  is  calculated 
using  the  schematic  in  the  same  figure.  The  differential  torque  created  by  the  force 
distributed  over  the  infinitesimal  length  d9R  is  dM  =  2FRd6Rsin{9).  This  can 
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Figure  50:  The  distributed  force  F  is  equivalent  to  the  moment  M.  R  is  the  radius 
of  the  ring  to  which  the  large  satellite  is  attached. 


be  readily  integrated  for  the  whole  ring  and  the  value  of  F  obtained. 

r  m 

M=  2 FR2  sin (0)dd  =>  F  =  —  (23) 

Jo 

where  the  moment  M  is  given  by  M  =  13 g  ■  m \  •  L  i ,  where  L\  is  known.  The  loads 
due  to  the  smaller  satellites  are  calculated  in  a  similar  manner. 

With  the  tools  presented  above  it  is  possible  to  study  any  PAD  design  encoded 
in  the  vector  x.  After  creating  the  PAD  structure  by  assembling  the  PAD’s  bare 
bones  and  the  several  Map  L-system  topologies,  the  tridimensional  PAD  model  is 
analyzed  via  the  FE  method  to  determine  the  fitness  of  the  individual.  For  this 
weight  minimization  problem  the  fitness  function  is  given  by: 

W(x)  +  Pa  +  Ps  (24) 


where  Pa  and  P$  are  the  penalization  terms  for  maximum  stress  and  maximum 
displacement,  respectively.  These  terms  are  given  by: 


p  —  5  femaX(amax-alirn,0)  _  ]^) 


Ps  =  5(' 


_,max((5ma£C  fiiiim-fi)  _ 


!)• 


(25) 


With  this  functional  form  for  the  penalization,  structures  for  which  the  maximum 
stress  and/or  maximum  displacement  are  marginally  higher  than  the  limiting  val¬ 
ues  are  not  very  penalized. 


52 


From  the  FE  method  solution  the  maximum  stress  and  the  maximum  dis¬ 
placement  can  be  determined.  The  stresses  on  the  Euler  beams  results  from  trac¬ 
tion/compression  and  also  from  bending.  To  avoid  complicated  calculations  to 
determine  the  actual  maximum  tension  on  each  beam,  a  conservative  value  for  the 
maximum  stress  in  the  structure  is  used: 


=  Ft 


max, bending]  T  tension] 


(26) 


where  the  values  amax,  bending,  and  a  tension  are  readily  available  in  the  COMSOL 
post-processing  functions. The  maximum  displacement  in  the  PAD  ( Smax )  is  also 
calculated  with  the  COMSOL  postprocessing  functions. 

The  mass  is  calculated  by  integrating  the  cross-sectional  area  along  the  beams 
length.  This  calculation  yields  a  conservative  value  for  the  PAD’s  mass  since  it 
takes  into  account  the  mass  in  the  intersecting  beams  twice,  as  exemplified  in 
figure  ??.  It  is  shown  ahead  this  can  cause  a  difference  of  more  than  15%  between 
the  mass  for  the  surrogate  model  and  the  mass  for  the  real  structure. 


Figure  51:  The  calculation  of  the  mass  is  conservative  because  the  mass  of  the 
cube  highlighted  in  the  intersection  of  the  two  beams  is  counted  twice. 


The  bioTOM  settings  for  the  evolutionary  algorithm  for  this  single  objective 
optimization  are  listed  below. 

•  There  are  100  individuals  per  population. 

•  The  initial  populations  were  generated  randomly. 

•  The  Genetic  Algorithm  ran  for  100  iterations  or  until  no  improvement  was 
observed. 

•  The  new  generations  were  formed  by  5%  elite  individuals,  80%  crossover 
offsprings  and  15%  mutated  individuals. 
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•  The  crossover,  selection  and  mutation  operators  are  as  explained  above. 

•  The  genome  is  encoded  in  vector  of  reals,  Y,  whose  elements  are  vary  between 
0  and  1,  just  as  it  was  described  above,  but  now  five  different  Map  L-system 
must  be  coded  instead  of  just  one  plus  the  thickness  for  the  PAD  bare  bones 
elements.  Figure  52  the  structure  of  Y. 


Map  L-system  1  Map  L-system  2  ^  Map  L-system  5  Bare  Bones  Thickness 


Figure  52:  The  genome  codification  for  the  PAD. 


The  convergence  history  for  the  best  design  is  depicted  in  figure  53. 

As  expected  the  curve  flattens  out  as  the  “best”  information  in  the  genetic 
pool  becomes  concentrated  in  one  individual.  After  70  generations  only  very  small 
improvements  are  observed  from  generation  to  generation. 

Figure  53  also  depicts  the  platforms  structures  for  the  best  individuals  in  in¬ 
termediate  generations.  An  important  conclusion  for  this  particular  case,  and 
perhaps  a  conclusion  that  goes  against  initial  conceptions,  is  that  the  best  struc¬ 
ture  is  a  rather  simple  one.  One  can  see  from  these  figures  a  constant  reduction  of 
the  number  of  beam  elements  for  the  best  structure  as  the  number  of  generations 
increases,  and  at  the  same  time  the  fine  tuning  of  the  beams  thickness.  Figure  54 
shows  the  beams  thicknesses  for  the  same  individuals  as  in  the  inserts  in  figure  53. 

The  optimal  structure  after  the  100  generations  is  depicted  in  figure  55.  The 
relevant  variables  for  this  individual  are  listed  below: 

•  the  mass,  masspAD  =  15.4  kg,  which  represents  a  reduction  of  30%  in  the 
mass  of  the  structure  when  compared  to  the  initial  design  (22  kg ); 

•  the  maximum  stress,  amax  =  294  MPa  is  below  the  300  MPa  as  required; 

•  and  finally  the  maximum  displacement,  Smax  =  1.6  cm,  is  also  below  the 
limiting  value  of  2  cm. 

In  order  to  validate  this  results  the  optimal  design  was  modeled  in  a  more 
realistic  manner.  The  beams  were  modeled  as  tridimensional  solid  elements  and 
the  satellites  were  included  in  the  model  to  simulate  the  loads  applied  to  the  PAD 
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Figure  53:  Evolution  of  the  best  PAD  design.  The  inserts  depict  the  structure  of 
the  platforms  for  some  selected  individuals,  all  of  them  the  best  of  their  respective 
generation. 


Figure  54:  The  thickness  of  the  beams  for  the  same  individuals  as  in  figure  53. 
The  color  represents  the  thickness  value  in  meters. 


more  accurately.  The  satellites  are  assumed  to  be  rigid  and  homogeneous  masses 
with  dimensions  such  that  their  centers  of  gravity  correspond  to  the  ones  from  real 
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Figure  55:  The  structure  of  the  3  platforms  for  the  best  individual  after  100 
generations,  and  the  assembled  geometry,  (a)  top  platform,  (b)  middle  platform, 
(c)  bottom  platform  and  (d)  assembled  structure. 


satellites.  The  simulated  geometry  is  depicted  in  figure  56. 

The  stresses  and  the  displacements  are  is  shown  in  figure  58.  In  this  case 
the  Von  Mises  stresses  were  used  as  the  failure  criterion.  The  results  show  some 
differences  with  respect  to  the  surrogate  model  but  in  general  the  agreement  in 
the  maximum  values  for  stress  and  displacement  and  the  location  in  the  structure 
where  they  occur  is  satisfactory,  as  seen  in  figures  58  and  57.  Furthermore,  the 
realistic  model  also  satisfies  all  the  displacement  and  stress  constraints.  As  for 
the  mass,  and  because  of  the  reasons  pointed  above,  its  value  decreased  to  13  kg, 
which  represents  a  mass  reduction  of  41%  from  the  initial  design. 

To  finalize  the  PAD  was  redesigned3  in  detail  following  what  was  learned  from 
the  bioTOM.  The  final  design  for  the  PAD  is  shown  in  figure  59.  In  the  end  the 
mass  of  the  structure  increased  slightly  to  13  kg  however  that  still  represents  a 
remarkable  decrease  of  41%  in  the  mass  with  respect  to  the  initial  design. 

§3.  MAIN  FINDINGS 

In  this  work,  a  novel  paradigm  for  aerospace  structural-topology  design  was 
introduced  and  assessed.  This  paradigm  is  inspired  on  the  cellular  division  of 
living  organisms,  which,  similarly  to  engineering  design,  depends  on  their  fitness 
to  survive  or  replicate.  The  following  is  a  summary  of  the  main  findings  of  this 
research: 


3The  author  thanks  again  to  Carol  Hude  for  providing  the  final  design. 
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Figure  56:  The  tridimensional  PAD  structure  with  the  satellites  assembled. 


1.  bioTOM  evolves  engineering  designs  that  are  higher-performing  than  the  designs 
obtained  with  Genetic  Algorithm.  In  test  cases  carried-out  comparing  the  pro¬ 
posed  bioTOM  and  Genetic  Algorithm,  bioTOM  consistently  evolved  higher¬ 
performing  designs.  The  main  reasons  for  this  improved  design  are  twofold: 
(a)  the  evolution-development  (EvoDevo)  paradigm,  and  (b)  bioTOM  gen¬ 
erates  a  much  larger  proportion  of  feasible  designs.  The  import  of  EvoDevo 
paradigm  derives  from  mainly  for  two  reasons,  it  provides  a  more  effective 
searching  of  the  design  space  and  it  provides  designs  with  better  resolution. 
The  former  advantage  can  be  easily  grasped  in  the  biological  context,  which 
inspired  bioTOM.  Indeed,  a  mutation  of  a  single  cell  would  provide  exceed¬ 
ingly  small  variation  in  the  design  to  be  an  effective  agent  of  change  in  the 
individual  fitness  (unless,  that  is,  the  design  is  so  unstable  that  a  vanish- 
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Figure  57:  The  finite  element  solution  simplified  (Euler  beams)  model  of  the  PAD. 
(a)  Stresses  [Pa]  calculated  using  equation  26,  (b)  the  deformed  structure  [m]. 


Figure  58:  The  finite  element  solution  for  the  tridimensional  model  of  the  PAD. 
(a)  The  Von  Mises  stresses  [PA],  (b)  the  deformed  structure  [m] 


ingly  small  variation  in  its  phenotype  would  lead  to  large  variation  in  the 
its  associated  fitness;  a  prospect  that  would  be  arguably  of  any  use  either 
in  biology  or  engineering.)  The  second  advantage  of  EvoDevo  allows  for  a 
more  effective  search  of  the  design  space  by  providing  fincreased  definition 
with  smaller  genes.  A  smaller  gene  in  its  turn  allows  for  a  more  efficient 
search  of  the  design  space.  The  other  main  advantage  of  bioTOM  over  Ge¬ 
netic  Algorithm  is  the  fact  that  a  much  larger  fraction  of  the  population  in 
each  generation  is  feasible  (test  case  1).  One  reason  for  that  stems  from  the 
network  structure  intrinsic  to  bioTOM  but  which  has  to  be  “learned”  by 
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Figure  59:  The  optimized  PAD  structure. 


trial-and-error  in  the  conventional  Genetic  Algorithm. 

2.  bioTOM  evolves  better  designs  than  Genetic  Algorithm  in  a  fraction  of  the  num¬ 
ber  of  generations  needed  by  Genetic  Algorithm.  This  advantage  follows  from 
the  more  effective  search  of  the  design  space  explained  above  using  bioTOM 
as  compared  with  the  Genetic  Algorithm.  Because  this  search  is  substan¬ 
tially  better  that  the  conventional  search  with  Genetic  Algorithm,  bioTOM 
can  obtain  higher  performing  designs  than  Genetic  Algorithm  in  one  order 
of  magnitude  less  generations  than  Genetic  Algorithm. 

3.  bioTOM  evolves  better  designs  than  Genetic  Algorithm  using  a  fraction  of  the 
number  of  genes  needed  by  Genetic  Algorithm.  Again,  this  is  a  consequence 
of  EvoDevo  model  used  in  bioTOM.  Indeed,  the  number  of  genes  in  living 
organisms  is  only  a  minute  fraction  of  the  number  of  cells  they  regulate. 
Thus  the  fact  that  the  DNA  encodes  a  program  that  develops  organisms, 
engenders  a  diversity  of  effective  solutions  that  can  be  found  and  evolved  in 
a  timely  manner.  This  advantage  was  shown  in  this  research  to  be  equally 
relevant  in  searching  for  optimal  engineering  designs  as  well. 
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4.  bioTOM  evolves  designs  that  are  as  high-performing  as  SIMP,  when  the  condi¬ 
tions  suit  SIMP,  but  provide  better  designs  than  SIMP  overall.  If  a  problem  is 
well  posed,  the  sensitivity  bounded  and  smooth,  the  fitness  landscape  convex 
and  the  optimal  topology  relatively  simple,  then  SIMP  will  converge  to  the 
optimal  design  as  the  mesh  is  refined.  The  problem  is  that  these  conditions 
are  seldom  met  in  real  applications — they  are  not  met  even  in  the  simplest  of 
the  test  cases  above.  It  is  nonetheless  reassuring  that  bioTOM  will  generate 
solutions  that  perform  as  well  as  SIMP  when  the  conditions  are  propitious 
for  SIMP  (test  case  2).  Overall,  bioTOM  will  perform  better  because  of 
the  increased  resolution  of  the  topology  as  explained  above  and  because  it 
searches  for  global  minima,  instead  of  local  minima  as  in  SIMP. 

5.  bioTOM  is  more  versatile  than  SIMP.  As  mentioned  in  the  previous  item,  SIMP 
requires  a  number  of  conditions  for  its  adequate  use.  Because  bioTOM  do 
not  require  these  same  conditions,  it  can  be  applied  to  a  more  varied  class  of 
problems  as  compared  to  SIMP.  For  instance,  bioTOM  does  not  require  the 
a  priori  definition  of  the  region,  called  design  domain,  over  which  the  design 
is  sought,  bioTOM  does  not  require  a  regular  fitness  function  or  constraint,  a 
maximum  stress  is  a  case  in  point,  bioTOM  does  not  simulate  empty  space, 
and  bioTOM  searches  for  global  minima,  instead  of  local  minima  in  SIMP. 
The  latter  is  particularly  relevant  for  engineering  applications,  where  the 
complexity  and  non-linearity  of  the  constraints  and  objective  functions  lend 
themselves  to  multi- minima  problems.  In  these  cases,  the  solution  using 
SIMP  would  greatly  depend  on  the  initial  condition. 

6.  bioTOM  is  more  easily  integrates  existing  COTS  software.  In  this  research, 
bioTOM  was  easily  integrated  with  three  different  software:  COMSOL,  NAS- 
TRAN  and  a  in-house  vortex  lattice  method  code.  By  contrast,  SIMP  re¬ 
quires  the  computation  of  sensitivities  and  these  computations  are  not  easily 
implemented  in  COTS  codes.  For  instance,  these  calculations  could  be  im¬ 
plemented  using  a  general  finite- differences  approach.  In  this  case,  the  price 
of  generality  would  be  poor  performance  at  best  and  inaccurate  and  diverg¬ 
ing  sequences  of  designs  at  worst.  Alternatively,  specific  and  exact  gradient 
calculation  could  be  sought.  However,  such  approach  would  be  hardly  inte- 
grable  with  existing  software,  as  it  would  require  programing  the  sensitivities 
for  each  objective  function  and  for  each  different  application.  That  could 
only  e  accomplished  with  access  to  the  source  code  of  the  COTS  software — a 
virtual  impossibility  in  most  cases. 
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7.  bioTOM  provides  better  solutions  for  engineering  design  than  conventional  bio- 
mimetics.  Instead  of  mimicking  extant  biological  structures  for  engineer¬ 
ing  design,  as  in  biomimetic  approaches,  the  proposed  methodology  seek  to 
model  and  simulate  the  EvoDevo  mechanisms  that  led  to  the  emergence  of 
those  structures  in  the  first  place.  This  distinction  is  critical  for  engineering 
design  because  the  requirements,  constraints  and  material  available  for  living 
organisms  are  utterly  distinct  from  engineering  design.  And  because  there 
are  some  applications  for  which  no  biological  model  exists,  spacecraft  design 
is  a  case  in  point. 

Although  more  studies  are  required  to  further  validate  and  refine  the  method¬ 
ology,  the  results  reported  in  this  work  clearly  demonstrate  that,  with  expected 
extensions  and  improvements,  the  methodology  can  play  a  key  role  in  the  design 
of  novel  and  better  performing  aerospacel  structures.  In  fact,  no  one  doubts  the 
prime  role  of  EvoDevo  in  the  evolution  of  the  species.  This  research  showed  that 
synthetic  EvoDevo  can  be  as  consequential  and  beneficial  for  engineering  design 
as  biological  EvoDevo  is  for  the  evolution  of  natural  organisms. 
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plan  to  pioneer  the  use  of  topology  optimization  in  the  development  of  light, 
high-performing  structures  for  the  satellite  bus  and  payload.. 
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4.  American  Society  for  Engineering  Education  Fellow  at  Wright-Patterson  Air 
Force  Research  Laboratory  2009  Air  Force  Summer  Faculty  Fellowship  Pro¬ 
gram.  M.  H.  Kobayashi  was  a  fellow  at  AFRL/WP  during  the  summer 
of  2009.  The  scientific  advisors  for  this  fellowship  were  Dr.  Raymond  M. 
Kolonay  and  Dr.  Gregory  W.  Reich,  both  at  AFLR/WP.  In  this  interaction, 
the  shape  and  topology  optimization  of  an  MAV  was  carried-out.  The  results 
of  this  interaction  will  be  presented  at  the  51st  AIAA/ASME/ASCE/AHS/ASC 
Structures,  Structural  Dynamics,  and  Materials  Conference  in  Orlando,  Florida. 
Graduate  student  A.T  LeBon  was  also  a  fellow  with  M.  H.  Kobayashi  at 
AFRL/WP  during  the  summer  of  2009.  During  this  fellowship,  Dr.  Kolonay 
and  M.  H.  Kobayashi  prepared  the  draft  of  a  successful  grant  with  AFOSR 
in  the  area  of  Structural  Mechanics.  This  interaction  also  provided  an  op¬ 
portunity  for  M.H.  Kobayashi  to  acquaint  himself  with  NASTRAN. 

5.  Air  Force  Research  Laboratory/Air  Vehicles  Directorate  (AFRL/RB)  2008 
Summer  Researcher  Program  Fellowship.  M.  H.  Kobayashi  was  a  fellow 
at  AFRL/WP  during  the  summer  of  2008.  In  this  interaction,  two  research 
projects  were  planned  and  executed.  The  first  project  aimed  at  optimizing 
the  mechanism  for  a  morphing  wing  technology.  The  mechanism  was  de¬ 
sign  and  the  results  were  published  as  an  AIAA  paper  for  the  50th  AIAA/ 
ASME/ASCE/AHS/ASC  Structures,  Structural  Dynamics,  and  Materials 
Conference,  Palm  Springs,  California.  The  second  project  investigated  the 
optimal  design  for  the  wing  box  of  a  generic  fighter.  The  results  of  this  re¬ 
search  were  published  in  The  Aeronautical  Journal  of  the  Royal  Aeronautical 
Society. 

6.  Class  on  Comparative  Biomechanics.  Another  important  interaction  during 
this  period  was  the  development  and  teaching  of  a  new  technical  elective 
course  on  Comparative  Biomechanics  for  students  of  both  Biology  and  En¬ 
gineering  Majors.  The  course  was  jointly  developed  by  M.  H.  Kobayashi 
and  M.  A.  Butler,  and  was  well  received  by  the  students,  with  fourteen  en¬ 
gineering  students  and  three  biology  students  enrolling  in  the  class — this 
enrollment  is  more  than  double  the  number  of  students  that  regularly  enroll 
in  technical  elective  classes  in  engineering. 

§7.  NEW  DISCOVERIES,  INVENTIONS  OR  PATENT  DISCLOSURES 

Provisional  patent  61/284,520  on  bioTOM. 
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§8.  HONORS/AWARDS 

1.  American  Society  for  Engineering  Education  Fellow  at  Wright-Patterson  Air 
Force  Research  Laboratory  2009  Air  Force  Summer  Faculty  Fellowship  Pro¬ 
gram. 

2.  Air  Force  Research  Laboratory/Air  Vehicles  Directorate  (AFRL/RB)  2008  Sum¬ 
mer  Researcher  Program  Fellowship. 
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